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Objectives 

The main objective of this research project was to estimate the thermal conductivity of 
cometary ices from computer simulations of model amorphous ices. This was divided into 
four specific tasks: 

1. Generating of samples of amorphous water ices at different microporosities; 

2. Comparing the resulting molecular structures of the ices with experimental results, for 
those densities where data was available; 

3. Calculating the thermal conductivities of liquid water and bulk amorphous ices and 
comparing these results with experimentally determined thermal conductivities; 

4. Investigating how the thermal conductivity of amorphous ice depends upon the micro- 
scopic porosity of the samples. 

Summary 

The thermal conductivity was found to be only weakly dependent on the microstructure of 
the amorphous ice. In general, the amorphous ices were found to have thermal conductivities 
of the same order of magnitude as liquid water. This is in contradiction to recent experi- 
mental estimates of the thermal conductivity of amorphous ice, and it is suggested that the 
extremely low value obtained experimentally is due to larger-scale defects in the ice, such as 
cracks, but is not an intrinsic property of the bulk amorphous ice. 
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Introduction 


It is thought that the ices in the region outside the solar system are representative of the 
protostellar" nebula from which the solar system was formed (Klinger, 1983). Comets and 
interstellar dust grains are thought to reflect the composition of the protostellar nebula, 
having existed in an unaltered state since the formation of the solar system. It is also 
thought that much of the water and biogenic material on the earth was delivered by cometary 
impacts after the initial formation of the planets. Thus, understanding the physical and 
chemical state of the icy interior of cometary bodies would provide not only a more detailed 
understanding of the early solar system, but would also yield information on the form and 
composition of biogenic elements as they were introduced to the earth (Chyba et ai, 1990). 

A great deal has been learned about astrophysical ices through infrared astronomical ob- 
servation of interstellar molecular clouds, planets, icy satellites and cometary bodies (Hagen 
et ai, 1979; Sandford and Allamandola, 1990b; Sandford and Allamandola, 1990a; Cruik- 
shank et al., 1993; Owen et al., 1993; Irvine et ai, 1996). Our understanding of astrophysical 
ices has been aided by ground-based experiments (Hallbrucker et ai, 1991; Klug et ai, 1991; 
Jenniskens and Blake, 1994; Devlin and Buch, 1995) and computational studies (Zhang and 
Buch, 1990a; Zhang and Buch, 1990b; Buch, 1990; Buch, 1992; Wilson et ai, 1995; Essmann 
and Geiger, 1995). Computational studies are particularly well suited for investigating con- 
ditions which cannot be obtained in the laboratory, such as the formation of cold (10-15 K) 
amorphous ice from cold water vapor. 

Cometary bodies which pass near the sun, however, are subject to solar radiation, giv- 
ing rise to possible structural changes in the icy mantle as well as out gassing of volatile 
constituents (Prialnik and Barnun, 1990; Notesco et ai, 1991; Notesco and Bar-Nun, 1997). 
The extent to which this processing occurs depends upon the thermal conductivity of the 
amorphous ice. In addition to its importance in cometary evolution, the thermal conductiv- 
ity is also important in the formation of ice on interstellar dust grains. Recent experiments 
have led to widely different values of the thermal conductivity, and the problem has not been 
resolved (Kouchi et ai, 1992). 


Methods 

In an MD simulation, the equations of motion of the system are solved numerically, and the 
properties of interest are obtained from averages over the resulting numerical trajectories. 
The majority of MD calculations are based on classical mechanical models, and the problem 
reduces to the solution of Newton’s equations of motion for a system of interacting atoms 
and/or molecules. A detailed descriptions of this method can be found in a number of basic 
texts (Allen and Tildesley, 1987; Hansen and McDonald, 1986). 

Several models are available to describe the interactions between water molecules. We 
have selected the TIP4P model (Jorgensen et ai, 1983) which has been shown to provide a 
good description of liquid water, as found through the comparison of the water-water pair 
correlation functions with the results of neutron scattering experiments (Jorgensen et ai, 
1983), as well as the liquid-vapor interface of water (Wilson et ai, 1987). The TIP4P model 
can also reproduce many physical properties of crystalline ices (Tse and Klein, 1988a), and 
the pressure disordered high and low density amorphous ices (Tse and Klein, 1988b). Related 
studies on solid amorphous ice phases have been done using the ST2 model of water (Poole 
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et al., 1992), and amorphous ice clusters using the TIPS2 model (Zhang and Buch, 1990a; 
Zhang and Buch, 1990b; Buch, 1990; Buch, 1992). 

The methodology for generating these structures has been described previously (Wilson 
et al. , 1995). Briefly, the TIP4P model of water was used throughout the simulations (Jor- 
gensen et al., 1983), which has been shown to reproduce correctly many of the structural 
properties of water and aqueous solutions. 343 water molecules were placed in a cubic box 
at a density of 0.94 g/cc for the Ida ice and 1.17 g/cc for the hda ice. The systems were 
simulated at a temperature of 300 K for 1.1 ns. Configurations were stored every 10 ps over 
the final 1.0 ns of the run. At 300 K, the correlation time for water is less than 10 ps, so this 
yields a set of 100 uncorrelated configurations at each density, for a total of 200 configura- 
tions. It is important to start with an uncorrelated set of configurations: At the temperatures 
of interest, 10-77 K, a trajectory started from a single configuration will undergo very few 
structural rearrangements on the time scale that can be followed in a computer simulation. 
Constructing averages over a set of uncorrelated configurations avoids this “initial condi- 
tions” problem. Each of the 200 configurations was slowly annealed to 10 K by cooling the 
configuration in 50 K increments over 50 ps followed by a constant temperature trajectory 
of 50 ps. 

Results 

Bulk Amorphous Ices 

We have carried out computer simulations of bulk amorphous ices at the densities of the two 
major polymorphs: 0.94 g/cc for low density amorphous (Ida) ice and 1.17 g/cc for the high 
density amorphous (hda) ice (Narten et al., 1976). The structure of the ices, as determined 
from the pair-correlation functions, was found to be in good agreement with the correlation 
functions obtained from neutron scattering experiments. This demonstrates that the water 
models used in our computer simulations provide a good description of the bulk amorphous 
ice phases. The neutron-weighted structure factors, G(r), of the low and high density ices 
computed from the bulk quenches are shown in Figure 1. These are in good agreement with 
experimentally determined structure factors (not shown) (Bellissent-Funel et al., 1987) The 
neutron-weighted pair correlation function, G(r), is related to the usual atom-atom radial 
distribution functions by G(r) = 4npr x (0.092Go_o(r)T0.422Go-i/( r ) + 0-486G//_//(r) — 1) 
where the numerical coefficients are scattering weights associated with the neutron scattering 
experiments. 

Detailed information about the nature of the hydrogen-bonding network of the ices can be 
obtained from the oxygen-oxygen radial distribution functions (RDF’s), shown in Figure 2. 
The oxygen-oxygen RDF of Ida ice is similar to that of liquid water, also shown in Figure 2. 
The location of the first and second peaks at 2.75 Aand 4.45 Ain the RDF are indicative of 
the tetrahedral ordering found in liquid water. This result is in good agreement with X-ray 
diffraction results on pressure induced Ida ice (not shown) (Bizid et al., 1987). The main 
difference between the computer simulation results and the X-ray data is that the first peak 
is much broader in the X-ray data, which suggests that the computer simulated ices may 
not be fully relaxed. 

In contrast, the first peak in the oxygen-oxygen RDF of hda ice is smaller and broader 
than for Ida ice, and the second peak is very broad and shifted to smaller distances, indi- 
cating a greatly distorted hydrogen bond network. The X-ray results on hda ice also show 
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Figure 1: Neutron- weighted scattering 

functions from computer simulations of Ida 
(solid line) and hda (dot-dashed line) ice. 



Figure 2: Oxygen-oxygen pair correlation 

functions for water at 300 K, Ida ice and hda 
ice. 


substantially less structure than the Ida ice. However, the agreement between the RDF of 
the computer simulated hda ice and the X-ray diffraction results (Bizid et a/., 1987) is only 
qualitative. In particular, the X-ray RDF of the hda ice exhibits two next-nearest neighbor 
peaks at 3.7 A and 4.65 A, of which the 3.7 A is due to the presence of water triplets with 
small 0-0-0 angles. This feature appears to be absent in simulated hda ice, although it 
is possible that the single broad peak observed here, which spans the region from 3.7 Ato 
4.7 A, is not completely relaxed, and would split into two peaks with additional annealing. 

These results demonstrate that the TIP4P model adequately represents the structural 
properties of bulk samples of the hda and Ida ices. The bulk Ida and hda samples where then 
used as the initial substrates in simulations to generate samples at different microporosities. 

Calculation of the Thermal Conductivity of Amorphous Ice 

The thermal conductivity, k, can be calculated from the long-time limit of the heat flux 
correlation function, 


* = 2 J^vt* JS2 5 < E mEMW) - *.(°)) 2 >. 

where ks is the Boltzmann constant, T and V are the temperature and volume of the system, 
respectively, and E,{t) and x x {t) are the energies and positions of molecule i at time t. The 
double sum is over all N molecules in the system, and < ... > denotes a statistical average. 

The heat flux correlation for liquid water is shown in Figure 3. The thermal conductivity 
of water is determined from the slopes of these curves in the long-time, liner regime. It is 
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found to be 0.19, 0.23 and 0.26 W-m -1 -K — 1 for temperatures of 290, 300 and 310 K, respec- 
tively. These values are between 2 and 3 times lower than the experimentally determined 
thermal conductivities of water at the same temperatures: 0.5917, 0.6096 and 0.6252 W-m -1 - 
K— 1. While this difference is fairly large, it is small compared to the differences between 
the measured thermal conductivity of water and experimental estimates of the thermal con- 
ductivity of amorphous ice. Thus, our computer simulations provide a sufficiently accurate 
description of the thermal transport properties of water that meaningful conclusions can be 
drawn from subsequent calculations of the thermal conductivity of the amorphous ice as a 

To create samples of microporous amor- 
phous ice, we carried out high temperature 
simulations at 300 K with fixed bulk densities 
of 0.9, 0.8 and 0.7 g/cc. Configurations from 
these calculations were quenched to 77 K to 
provided a set of initial conditions for calcu- 
lations of the thermal conductivity. To inves- 
tigate how different microporous structures 
might affect the thermal conductivity, we de- 
vised two methods for generating such struc- 
tures. In the first method, we took a set of Ida 
ice configurations and randomly deleted wa- 
ter molecules until the density was 0.7 g/cc. 
The second procedure was to delete water 
molecules around a fixed center until a den- 
sity of 0.7 g/cc was reached. This allows us 
to investigate the possible extremes of the microporous structures: random deletion of water 
molecules results in a fairly uniform distribution of small cavities in the structure whereas 
deleting around a fixed center yields much larger cavities. 

The heat flux correlation function for the 
two microporous samples of amorphous ice 
are shown in Figure 4. These are qualitatively 
similar to the results for liquid water, shown 
if Figure 3, in particular, both the liquid and 
amorphous samples reach the linear regime 
on the same time scales and the curves have 
slopes that are of the same order of magni- 
tude. The main difference between the liquid 
and amorphous figures is that the heat flux 
correlation of for the amorphous substrates 
makes a more abrupt transition to the linear 
regime. This is due to the lack of relaxation 
in the energy flux via mass transport in the 
amorphous ice samples over the time scales 
of the computer simulation. What is most 
significant is that the results of random deletion of water molecules, giving rise to a large 
number of smaller micropores, is not qualitatively different than the results for deletion of 



Figure 4: Heat flux current correlation for mi- 
croporous amorphous ices at 0.7 g/cc. 


function of the degree of its microporosity. 
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Figure 3: Heat flux current correlation function 
for bulk liquid water at 290 K, 300 K and 310 K. 
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Thermal Conductivities of Water and Amorphous Ices 



T(K) 

MD 

K (W- 

exptl 

-m _ 1 -K _1 ) 

Water 

290 

0.19±0.02 

0.5917“ 


300 

0.23±0.02 

0.6096“ 


310 

0.26±0.01 

0.6252“ 

Ice Ih 6 

100 

1.49 

5.94 


260 

0.94 

2.35 

LDA 

93 

0.15 

4xl0 -5 (80-100 K) c 

HDA 

93 

0.35 


microporous^ 




0.9 g/cc 

86 

0.12±0.02 


0.8 g/cc 

85- 

0.10±0.02 


0.7 g/cc 

86 

0.10, 0.23' 

(87 K), 0.14' (92 K) 


“Ramires,et al., J. Phys. Chem. Ref. Data (1995) 24, 1377-1381. 

6 Ice Ih results are from: Inoue, et al., J. Chem. Phys. (1996) 102, 9569-9577. 

c Kouchi, et al., Astrophysical J. 388, L73-L76. 

d micropores created by annealing at fixed density from liquid. 

e micropores created by random deletion of water molecules followed by quenching. 

' micropores created by deletion around a fixed center followed by quenching. 


water molecules around a single center, creating a single, large void in the sample. Conse- 
quently, even at the low density of 0.7 g/cc, we find no evidence for a greatly reduced thermal 
conductivity in either extreme of a few, large micropores or uniformly dispersed micropores. 

The computed thermal conductivities of the various water and amorphous ice runs simu- 
lations are summarized in Table 1, as well as the computed thermal conductivities of the hda 
and Ida ice. They are of the same order of magnitude as the results for liquid water. This 
result is reasonable because the structure of amorphous ice is quite similar to liquid water. 
The thermal conductivity of crystalline ice Ih at 100 K is 5.95 W-m -1 -K -1 , over an order of 
magnitude larger than the thermal conductivity of amorphous ice at similar temperatures. 

The main result is that the thermal conductivity of the microporous ices appear to be 
of the same order of magnitude as the bulk amorphous ices. The thermal conductivity does 
not decrease dramatically as the density of the amorphous ice is reduced from 0.94 g/cc 
to 0.7 g/cc. Since the thermal conductivities for ice samples containing large cavities and 
uniformly distributed small cavities differ by only a factor of 2 we conclude that the bulk 
thermal conductivity is not very sensitive to the distribution of the micropores. Therefore, 
the microporosity of the ice does not explain the extremely low value of thermal conductivity 
that has been reported for vapor deposited amorphous ice. 


Significance 

Since all hypotheses about the role of comets in the origin of life and chemical evolution of 
the solar system make explicit or implicit assumptions about the thermal conductivity of 
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the cometary ice, the results of our calculations provide crucial information for the critical 
evaluation of these hypotheses. Recent experiments on laboratory analogs of cometary ice 
have yielded values of the thermal conductivity that are several orders of magnitude lower 
than the thermal conductivity of bulk amorphous ices. If this is true, then it implies that 
the interiors of comets are not heated at the perihelion of their orbits. Consequently, the 
interiors of comets might not have undergone any significant amount of thermal processing 
and their composition might reflect the composition of the protostellar nebula. 

Whereas the laboratory analogs of cometary ice are prepared by slowly depositing water 
vapor onto a cold substrate, the bulk amorphous ices are prepared by applying extreme 
pressure to crystalline ice samples until the underlying lattice is disrupted. It was originally 
thought that the low measured value of the thermal conductivity might be due to the micro- 
scopic porosity of the ice sample. As the computed thermal conductivities of several samples 
of microporous ice do not appear to be extremely sensitive to the microscale structure of the 
ice, we conclude that an extremely low thermal conductivity is not a property of amorphous 
ice, even at densities as low as 0.7 g/cc considering several extremes of how the microscale 
structure of the ice is induced. Rather, our results indicate that the thermal conductivity of 
microporous, amorphous ice samples are similar to bulk amorphous ices. Consequently, it is 
likely that the thermal conductivity of astrophysical ice is of the same order of magnitude 
as normal ice. 

It is possible that the extremely low thermal conductivity measured in the vapor deposi- 
tion experiments is due to larger scale defects, such as cracks. In this case, the macroscopic 
thermal conductivity is determined by the heat flux across such defects and not the bulk 
thermal conductivity of the bulk amorphous ice. In the laboratory, the initial density of the 
ice will depend upon a number of factors including the deposition rates and the temperature 
of the substrate (Jenniskens and Blake, 1994), and amorphous ice layers from smooth and 
dense (Jenniskens and Blake, 1994) to extremely porous and “fluffy” (Laufer et ai, 1987) 
can be obtained with different initial conditions. 

Publications: “The Thermal Conductivity of Cometary Ice”, M. A. Wilson, and A. Poho- 
rille, (to be published, 1998). 

Keywords: thermal conductivity, microporosity, amorphous ice, cometary ices 
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